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Abstract 

We consider resequencing studies of associated loci and the problem of prioritizing se¬ 
quence variants for functional follow-up. Working within the multivariate linear regression 
framework helps us to account for the joint effects of multiple genes; and adopting a Bayesian 
approach leads to posterior probabilities that coherently incorporate all information about the 
variants' function. We describe two novel prior distributions that facilitate learning the role of 
each variable site by borrowing evidence across phenotypes and across mutations in the same 
gene. We illustrate their potential advantages with simulations and re-analyzing a dataset of 
sequencing variants. 


1 Introduction 


Genome-Wide Association Studies (GWAS) have allowed human geneticists to compile a rather 
long list of loci where DNA variation appears to be reproducibly associated to phenotypic vari¬ 
ability (|NHGRI 20151. While these might represent only a subset of the portion of the genome 
that is important for the traits under study ( Manolio et a?.||2009| l, there is little doubt that under¬ 
standing the characteristics and mechanisms of functional variants at these loci is a necessary 
next step. As resequencing becomes ever more affordable, follow-up investigations of GWAS loci 
often start with a comprehensive catalogue of their genetic variants in a sample of thousands of 
individuals, raising the question of how to sort through these results. 

Among the many challenges, let us discuss two. First, common variants are often correlated 
and it is difficult to distinguish their roles without accounting for the broader genetic background 
of the individuals who carry them. Second, rare variants are present in a small enough portion 
of the sample that statistical statements become impossible. With this in mind, it has been noted 
that (a) it is important to account for correlation between variants to obtain useful ranking; 
(b) we should increasingly be able to take advantage of the information gathered through other 
studies; and (c) Bayesian models provide a principled approach to guide variant prioritization. 
To adequately select among variants in the same locus (defined as a genomic region that might 
encompass multiple genes but that corresponds to the same association signal in a GWAS study), 
researchers have resorted to model selection approaches (Valdar et fl/.|2012| or approximations of 
the joint distribution of univariate test statistics ( [Faye et fl/.||2013 |Hormozdiari et ^|2014[). Prior 
information on variant annotation has been incorporated in models for eQTL (Veyrieras et al. 
2008 1 and more recently for general traits (Pickrell 2014| Kichaev et al. 2014| [Ghung et al. 2014 1 , 
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and annotation programs increasingly attempt to include information on identified genetic loci 
(|Wang et al. 2010[|. Prioritization often relies on Bayes's theorem, and Bayesian methods have 


received renewed attention in the context of GWAS data analysis (Guan and Stephens [2011 


Peltola et flZ.|2012a|b|l, genomic prediction (|Gianola|2013|l, and the evaluation of heritability (IZhou 

effll.||2013| 


In this context, we explore the advantages of a careful specification of the prior distributions 
on variants, by allowing sharing of information across multiple phenotypes and across neighbor¬ 
ing rare variants. We are motivated by the analysis of an exome resequencing study ([Service et al. 


2014|| in which individual level data is available for exomic variants at multiple genomic loci that 


have demonstrated evidence in GWAS of association to lipid traits. By design, the vast majority 
of measured variants is coding or in UTR, that is, in portions of the genome with high prior 
probability of harboring functional mutations. Annotation can help distinguish the role of syn¬ 
onymous variants, and conservation scores can be used to predict the effect of nonsynonymous 
ones; but annotation cannot be used to discount the importance of a large number of non-coding 
variants that one can expect to occur in a whole genome sequencing dataset. Measures on levels 
of High Density Lipoprotein (HDL), Low Density Lipoprotein (LDL), and triglycerides (TG) are 
available for the study subjects, and we are interested in capitalizing on the multidimensional 
nature of the phenotype. Prior analyses of this dataset ( Service et fl/.||2014| Bogdan et al.\m pressj l 
have illustrated the importance and the challenges of multivariate linear models, and we ex¬ 
plore here the advantages offered by carefully selecting priors for Bayesian models. Abstracting 
from the specifics of this dataset, we show how hierarchical prior distributions can be adapted 
to learn about the functionality of a variant by (i) looking across multiple phenotypes and (ii) 
aggregating the effects of multiple rare variants in the same gene. Since the power of Bayesian 
methods in borrowing information is well known, it is not surprising that others have explored 
their application in this context. For example, Yi et fl/.||(2011|| ill ustrate the use of priors to model 
a group effect for multiple rare variants, while Stephens ( [2013 1 describes models for the analysis 
of multiple traits. Our approach, however, is distinct from others in that it strives to achieve all 
of the following: (1) constructing a multivariate linear model that simultaneously accounts for 
the contributions of multiple genes and genomic loci; (2) providing inference on variant-specific 
effects—^while linking information across traits and genomic sites; and (3) accounting for the 
large number of variants tested, effectively enacting a form of multiple comparison adjustment. 

This paper is organized as follows. We devote Section to introducing the novel priors 
in the context of the genetic model, using an approximation of the posterior distribution to 
illustrate their inferential implications. Section [^describes the MGMG scheme used to sample the 
posterior, the setting used for simulations, and the criteria for comparison of methods. Section 
presents the results of simulation studies highlighting the potential of our proposal, as well as 
the description of the analysis of the motivating dataset. 


2 Prior distributions on genetic variants 

One characteristic of a genetic study based on resequencmg, as contrasted to genotypmg, is that 
researchers aim to collect a comprehensive catalogue of all genetic variants. This has implications 
for the statistical models used to analyze the data and the prior assumptions. Let n be the number 
of subjects in the study and p the number of polymorphic sites assayed. We will use y, to indicate 
the phenotypic value of subject i and the genotype of this subject at variant v (typically coded 
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as minor allele count). The simplest genetic model for a heritable phenotype is of the form 

Vi = 

keT 

where incapsulate all nongenetic effects and for k ^ T represent the contributions of a 
set T of genes that act additively and independently. Without loss of generality and following 
a standard practice in GWAS, we will assume that the effects of nongenetic determinants of the 
phenotypes have been regressed out from y, so that can be considered independent 'error' 
terms. Let us assume that the genetic effects are a linear function of minor allele counts so that 

yi — ^ f^v^iv G (1) 

veV 

for a set V of causal variants with e, iid Af{0,l/p)- Although this assumption is substantial, it 
only has the role of simplifying notation. While 0 represents the true genetic architecture of the 
trait, the membership of V is unknown in a typical association study, so the relation between the 
phenotype and genetic variants is expressed as 

Vi = E G ~ AT f 0, ^ y (2) 

U = 1 \ P J 


summing over all variable sites and with the understanding that only an (unknown) subset 
of ^ = {fii,..., fip) is different from 0. Below we will use the compact matrix notation y = 
Xj8 + e. Using (|^ to describe the relation between traits and genotypes depends heavily on the 
assumption that a resequencing study assays all variants. In GWAS, on the other hand, causal 
variants might be untyped, which means their contributions are partially captured by correlated 
variants and partially included in the error term. It would still be meaningful in that context 
to use a linear model to link phenotype and genotypes. However, in GWAS, the errors cannot 
be assumed independent, and the interpretation of the coefficients of X—as well as their prior 
distribution—is substantially more complicated. We note that mixed effects models can be used 
to address the first concern ( Kang et flZ.||2010[ |. 

The parameters in model (|^ are ^ and p; we now focus on their prior distribution. Follow¬ 
ing standard practice, we take p ~ Gamma(ap,Ap). (See Guan and Stephens (2011|l for another 
approach that specifically targets GWAS and relies on heritability information.) On the vector 
we want a prior that reflects our model selection goals and our understanding of the genetic 
architecture. There are several aspects to consider: (a) given the exhaustive nature of the geno- 
typrng process, we believe that most of the variants available do not directly influence the trait; 
(b) it seems reasonable that a variant that influences one trait (so that its effect size is definitely 
not zero) might also influence other traits; and finally (c) it appears likely that if a rare variant 
influences the outcome, other nearby rare variants might also have an effect. Our main goal is to 
describe prior distributions on that incorporate these beliefs. We start by recalling one class of 
priors that reflect (a) and then move on to generalizations that account for the sharing of infor¬ 
mation implied by (b) and (c). In what follows, we assume that the allele counts in the column 
of X have been standardized to have mean zero and variance one. 


2.1 Priors incorporating sparsity 

The prior belief that only a fraction of the typed variants has an effect on the phenotype is 
but one instance of what is a common assumption in high-dimensional statistics, i.e. that the 
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parameter of interest is sparse. To specify a prior on that gives positive probability to 


vectors with a number of coordinates equal to zero, we rely on a construction by George and 


McCulloch ( 1993| and introduce a vector of indicator variables Z such that = 0 implies /Sj, = 0 
The Zd are iid Bernoulli with parameter co, which governs the sparsity of the model and has a 
prior. Let indicate the collection of elements of )S corresponding to nonzero 
elements of Z, and let Xz be the corresponding columns of X. It has been found useful to assume 

()8z|Z,p,t) ~ Af (o, , where Ez is a known matrix and t ~ Unif(Ti, T 2 ) links the error 

variance to the size of the coefficients. In the literature, Ez mainly has one of two forms: 
J|z| (the identity matrix of size |Z|, where |Z| indicates the number of nonzero components of 


the vector Z) or n{X^Xz)^^, which is referred to as the g-prior (Zellner 19861 (and is a viable 
choice only when |Z| < n). Various views on the choice of Ez have been put forth (jChipman 
et fl?.||2001[ Heaton and Scott]|2010| Guan and Stephen^|2011[ |, but the strongest argument for the 
g-prior is that it provides computational benefits (see below). For either choice of Ez, all of its 
diagonal entries are equal, resulting in an equal prior variance for each of the /S^. Given the 
standardization of the columns of X, this implies that the original effect sizes are expected to be 
larger for rare variants than for common variants, which is reasonable. 

Sparse prior _ 


Bern(a;) 


cu~Beta(Aa;,Ba;) 


jO Gamma(ap, An) \p 


h 


T T" Unif[Ti,T2] 

)3z~A^(0,^Ez) 




Figure 1: Schematic representation of the Sparse prior distribution on j8. Hyperparameters are 
indicated in blue. The red portion describes all random objects and their dependency structure; 
unless explicitly indicated, random variables are independent. Boxes identify variables that share 
one of the distributions depicted in black. 

One of the advantages of the prior summarized in Figure Q is that the derived posterior 
distribution can be analytically integrated with respect to co, and p. While a MCMG is still 
needed to fully explore the posterior and carry out inference, we can rely on a collapsed Gibbs 
sampler that focuses only on t and the indicator variables Z. This reduces the computation at 
each iteration and improves its convergence rate ( |Liu 1994[ |. The prior densities for t and Z 
will be denoted, respectively, as /t(t) and /z(Z) —the latter being easily obtained from the beta- 
binomial distribution assumed for | Z |. As shown in the appendix, integrating jS and p out gives 
the marginal posterior density 


/z,T(Z,T|y) oc/^(t)/z(Z)(Ap + 



(f+v) det(Qz)^^^ 

Tlzidet(Ez)i/2' 


(3) 


where = X^Xz -|- and S| = y^y — y^XzOzX^y. Ghoosing Ez as in the g-prior leads 

to a simplification of the ratio in (|^, thereby avoiding the evaluation of one determinant at each 
iteration. 

Despite the need to evaluate numerically interesting summaries of the posterior of Z, we 
obtained an approximation (whose derivation and applicability is described in the appendix) to 
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gain a general understanding of how hyperparameters and data contribute to the final inferen¬ 
tial results. Specifically, we focus on the posterior expected value of Z^,, indicator of variant v, 
conditional on the indicators of all other variants Z[_j,]. In the case of orthogonal regressors, this 
expectation can be approximated as 


E[Zj;|Zr „,,T,y]^ 


T^Jn 


.-B,. 


p- |Z[^]| -1 


A. 


d-d 


(1 - vl) 


n/2 


(4) 


where pi, = xjy/ \/nj^ is approximately the correlation between variant v and the trait. From 
Q, one gathers that increasing |Z[_j,]|, which is the number of variants already used to explain 
the trait, increases the chance of an additional variant v to be considered relevant. This is a 
consequence of the fact that the parameter co, which describes the sparsity of j8 and hence the 
degree of poligenicity of the trait, is learned from the data (rather than set at a predetermined 
value). When a large number of variants have been found relevant, the trait is estimated to be 
highly poligenic and hence it is judged more likely that an additional variant might contribute to 
its variability. On the other hand, augmenting the total number of genotyped sites p will make it 
harder for any specific variant v to be judged important; this is adjusting for the look-ever 3 rwhere 
effect, an important step in gene mapping studies. 

Now that we have introduced this basic framework, we can consider modifications that facil¬ 
itate learning about the role of a variant across multiple traits and in the context of neighboring 
sites. We start with the first problem. 


2.2 Learning across traits 


One of the characteristics of current genetic datasets is the increased availability of multidimen¬ 
sional phenotypes. This is due partly to the automation with which many traits are measured 
and partly to the increased awareness that precise phenotypic measurements are needed to make 
progress in our understanding of the underlying biological pathways. Having records of mul¬ 
tiple traits in the same dataset allows for cross-pollination of genetic information. On the one 
hand, if a genetic variant is functional, it can be expected to impact more than one phenotype. 
On the other hand, even if noise in one phenotype makes it hard to distinguish the predictive 
power of a causal variant from that of a non-causal neighboring variant, it is much less likely 
that multiple traits would have noise correlated in such a way that causal and non-causal variants 
are indistinguishable for all of them. With this in mind, let us generalize the variant selection 
problem described in the previous section to handle multiple traits. 

Extending the notation, let jt be the standardized values for trait t, the coefficients of 
X in the mean E[yt], and Zf the corresponding indicator vector. We organize these by column 
in a n X q matrix Y, a p x q matrix jS, and a p x q matrix Z. Also let denote the 

number of traits associated with variant v, let be the entries of corresponding to entries 
equal to one in Zf, and let Xz, be the corresponding columns of X. The data-generating model 

is jt = ^Zt^z, + with et ~ A/^ ^0, and the priors on pt are simple extensions of the one 

used previously: pt iid Gamma(ap,Ap). Note that this model assumes that, conditionally on the 
genetic variants that influence them, the traits are independent; specifically, there are no shared 
environmental effects. This assumption might or might not be appropriate depending on context, 
but the prior distribution on jS that we are about to describe can be used also for models that do 
not rely on this assumption. 

We want a prior for )8 that continues to enforce sparsity but that allows learning about the 
role of a variant across traits. One possibility, first proposed by Jia and Xu (2007[|, is to introduce 
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a variant-specific probability of functional effect constant across traits and a priori indepen¬ 
dent with i/d ~ Beta(Ai„ B^;), where A-u and B^ can capture annotation information. Following 
the setup of the previous section, we then take (Zt,f|vi,) independent Bernoulli (v^,), set ^vt = 0 
whenever Z^t = 0, and let t) be independent across t with distribution Af{0, ^EzJ. 

As before, t Unif(Ti,T 2 ). 

As detailed in the appendix, we can derive an approximation analogous to Q: 


E[Z,t|Z[_(,,)],T, Y]-i « 1 + ^ (1 - nit) 


2 ^^/2 


(5) 


where t]vt = 


x?yt 

V^yJyt 


and Sj, [„f] = tallies the number of phenotypes for which the variant 


V has been judged relevant. This highlights a consequence of the selected prior distribution: as 
the total number of phenotypes q here has taken the role of p in Q, the role of each variant 
is judged not in reference to all the other variants but only in comparison to the effect of the 
same variant across traits. In other words, while there is learning across phenotypes, there is 
no adjustment for the multiplicity of queried variants. Bottolo et al. ( 2011[ | previously observed 
that sparsity of Zf could not be controlled by specification of the priors in this approach, and 
proposed letting Zyt have Bernoulli parameter VyOOt with independent priors on each factor. 

We propose a different remedy by introducing another layer in the hierarchical priors. 
Let W be a vector of indicator variables of length p: if W-o = 0, then Vi, = 0; if W^, = 1, 
v-v ~ Beta(Ai„Bz,). We take V\f-v iid Bernoulli(cuw) with ccw ~ Beta(Aw,B^y); the (Z^fjvz,) are 
independent Bernoulli(vy), as before. The schematic in Figure [^summarizes this prior proposal. 


Across Traits prior 



Figure 2: Schematic representation of the Across Traits prior distribution on jS. Hyperparameters 
are indicated in blue. The red portion describes all random objects and their dependency struc¬ 
ture; unless explicitly indicated, random variables are independent. Boxes identify variables that 
share one of the distributions depicted in black. 

The existence of a specific parameter for each site v allows variation in the average number 
of impacted traits per variant; some variants can be highly pleiotropic, while others are relevant 
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for one trait only. The sparsity parameter coy^j is once again estimated from the data, allowing 
for multiplicity adjustment. The introduction of W effectively specifies a hierarchical prior on 
Vp-, among the many possible ways to accomplish this, the one we adopt emphasizes the 
role of the sparsity parameter ooyj and is easily mterpretable. The appendix presents an indicative 
approximation of the posterior conditional expected values of Wy similar to Q and (|^. It de¬ 
pends on all phenotypes (enabling learning across traits), but the total number of variants p has 
again become the leading factor for effective multiplicity correction. We will refer to this prior 
as learning Across Traits. We include the first proposal in some comparison studies, indicating 
it as the Unadjusted approach to emphasize the fact that it does not include an adjustment for 
multiplicity. 

This may be an appropriate point at which to clarify the relation between the prior we 
are proposing and the traditional investigation of pleiotropy versus coincident linkage. The 
latter terminology derives from linkage studies, where the nature of the signal is such that the 
localization of variants with possible pleiotropic effects is possible only up to a certain genomic 
interval. This interval might contain one variant affecting multiple traits or contain different 
variants, each affecting a subgroup of the traits. First, it is worth noting that in this paper we 
are working in the context of association studies, which allow for a much finer resolution than 
linkage studies. The occurrence of multiple variants affecting multiple traits within the same 
linkage disequilibrium block is less likely, given that LD blocks are shorter than linkage regions. 
Secondly, ours is a fixed effects model using sequence data. We are aiming to estimate the specific 
effect of each variant rather than simply identifying a locus with a random effects model. Our 
framework, then, automatically considers two options: one variant affecting multiple traits or 
multiple variants affecting separate traits. The choice between these two alternatives is made 
on the basis of the posterior probability of the two models. This being said, it is important to 
recall that if two neighboring variants in LD affect two separate traits, the posterior probabilities 
of the two alternative models might be similar. The prior we introduce favors pleiotropy in the 
sense that it recognizes as likely that some variants affect multiple genes, but it does not exclude 
the alternative explanation, allowing the data to tilt the posterior in either direction. We have 
investigated this with simulations in the supplementary material. 

2.3 Learning across sites 

We now consider another form of 'learning from the experience of others' to improve our ability 
to identify functional variants. We focus on rare variants, which are observed in just a handful 
of subjects and for which it might be impossible to estimate individual effects. It is reasonable 
to assume that if one rare variant in a gene has an impact on a trait, other rare variants in the 
same gene also might be functional; with an appropriate hierarchical prior we might increase 
our ability to learn the effects of these variants. Of course a similar assumption might also be 
reasonable for common variants, but given that we observe these in a sufficiently large sample, 
we aim to estimate their individual effect without convolving their signal with that of others. 

The data-generating model is again (|^. We define r groups of variants, and we use 7(u) 
to indicate the group to which variant v belongs. Let G = (Gi,..., Gr) be a vector of indica¬ 
tor variables associated to the groups; we use these to link information from different variants. 
Specifically, if Gg = 0, then the proportion Vg of causal variants in group g is equal to zero; 
otherwise, Vg ~ Beta{Ag,Bg) (setting Vg = 1 for groups comprised of only one variant). The 
variant-specific indicators are iid Bernoulli with parameter Similarly to prior specifi¬ 

cations, (G^Iccg) are iid Bernoulli(6i;G) with coc ~ BetajAc,B q). This results in the partially 
exchangeable prior on jS represented in Figure the parameter Vg allows sharing information on 
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functionality across all variants in the same group. 


Across Sites /3 prior 



Figure 3: Schematic representation of the Across Sites prior distribution on j8. Hyperparameters 
are indicated in blue. The red portion describes all random objects and their dependency struc¬ 
ture; unless explicitly indicated, random variables are independent. Boxes identify variables that 
share one of the distributions depicted in black. 


As described in the appendix, the posterior conditional probability that a variant v belongs to 
the model depends on the overall number of groups, the number of groups considered relevant, 
and the number Sg = of variants in the same group that are deemed functional. The 

prior distribution in Figure 1^ which we refer to as learning Across Sites, allows one to achieve 
an effect similar to that of burden tests, while still providing some variant-specific information 
(which is in contrast, for example, to the proposal in Yi et al. (|2011|>). 


3 Methods 

3.1 MCMC sampling 

While we have resorted to some analytical approximation for expository convenience, we explore 
the posterior distribution with MCMC. As previously mentioned, we can focus on sampling t 
and all indicator variables. We use a Metropolis-within-Gibbs scheme, with the proposal distri¬ 
butions described below. For t, the common practice of using a truncated Gaussian works well. 
The discrete indicator variables pose a greater challenge, even though having integrated out jS 
allows us to work with sample space of fixed dimension, eliminating the need for a reversible 
jump MGMG. When there is only one layer of indicator variables Z, the proposal consists of first 
choosing with equal probability whether to add or remove a variant and then choosing uniformly 
among the candidate variants the one for which to propose a change of status. If the prior dis¬ 
tribution is described using higher level indicators as well, then proposed changes to both levels 
must be consistent. If an entry of W is changed from one to zero, the associated entries of Z also 
have to be zeroed; when proposing to change an entry of W from zero to one, the associated 
entries of Z are selected from the prior marginal. Additionally, there are proposal moves that 
leave W unchanged but then randomly select one of its nonzero entries and draw a proposal for 
the associated entries of Z in a fashion analogous to that described previously. Details of the 
algorithm are in the supplementary material. 

These simple proposal distributions will have trouble in two situations. The most common 
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is when two or more variants are strongly associated with a phenotype but are also strongly 
correlated with each other due to LD. Any specific Markov chain will tend to include one of the 
variants in the model, leaving out the rest. Another problematic situation is when the effects 
of two variants on a phenotype depend upon each other, so neither variant is likely to enter 
the model by itself, even if their joint inclusion would be favored by the posterior distribution. 
Others ( |Guan and Stephens|2011 Peltola et fl?.|2012a]b] l have described proposal distributions that 
overcome these difficulties and that can be reasonably applied to our setting—even though we 
do not investigate this in detail, focusing on the description of novel priors. 

The average Z of realized values of Z can be used to summarize the evidence in favor of each 
variant. Given its practical importance, the basic convergence checks incorporated in our package 
are based on Z. By default, the R code distributed in the package ptycho starts four chains 
from different points, runs each chain for a specified number of MGMG iterations, computes the 
averages for each chain separately, and then checks the range AZ of these averages. Details on 
the MGMG can be found in the supplementary material. The algorithm is implemented in the R 
package ptycho ( |Stell|2015^ . 


3.2 Evaluation of variant selection performance 

To investigate the performance of the proposed priors, we apply them to simulated and real data. 
The posterior distribution can be summarized in multiple ways. One can look for the indicator 
configuration that receives the highest posterior, for example, or make marginal inference on 
each variant. Both computational and robustness considerations make it practical to rely on 
posterior averages Z^t for comparisons. In the Bayesian models, then, we consider selecting 
variant v for trait t if the posterior average Z^t is larger than a certain threshold ^ e (0,1): 
iSf = {u : Zyf > ^}. For benchmarking purposes, we will also analyze the datasets with some 
non-Bayesian approaches. Specifically we will consider (a) the Lasso (Tibshirani 1996(b) a 
set of univariate linear regressions (one for each trait and variant), leading to f-statistics used 
to test the hypotheses of no association H^t ■ f^vt = 0 with multiplicity adjustment for the pq 


hypotheses via the Benjamini-Hochberg (BH) procedure at level a (Benjamin! and Hochberg 
1995| l; and (c) multivariate regression including all possible variants, with subsequent tests on 
the pq null hypotheses for each coefficient incorporating adjustment via the BH procedure at 
level a. The set of selected variants is equivalent in (a) to the set of estimated nonzero coefficients 
and in (b) and (c) to the set of variants for which the Hut : fivt = 0 are rejected. We will refer to 
these approaches as (a) Lasso, (b) BH marginal, and (c) BH full. 

The threshold ^ for Bayesian selection, the penalty of the Lasso, and the level a of BH can 
all be considered tuning parameters. We will compare the results of different procedures as 
these are varied (see details in the supplementary material). We base our comparison on an 
empirical evaluation of power and FDR associated with the different methods. Specifically, for 
each simulation and each method of analysis, we calculate the proportion of causal variants 
that are identified and the proportion of selected variants that are in fact false discoveries. The 
average of these values across multiple simulations is what we refer to as power and FDR in 
the results. The Bayesian methods also provide an estimate of FDR: if Z^t is approximately the 
probability that variant v is causal for trait t, then the mean of (1 — Z^t) over the selected variants 
is the Bayesian False Discovery Rate. We let BFDR denote this mean and explore how well it 
approximates (or not) the realized FDP, evaluated across all traits and variants. 
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3.3 Genotype and phenotype data 


Our work has been partially motivated by a resequencing study: Service et al. ( 2014[ | analyzed tar¬ 
geted exome resequencmg data for 17 loci in subjects of Finnish descent (from the 1966 Northern 
Finland Birth Cohort (NFBC) and the Finland-United States Investigation of NIDDM Genetics 
study (FUSION)). While the original study considered six quantitative metabolic traits, we focus 
here on the fasting levels of High Density Lipoprotein (HDL), Low Density Lipoprotein (LDL), 
and triglycerides (TG), transformed and adjusted for confounders as in the initial analyses (see 
supplementary material). The genotype data was obtained by sequencing the coding regions of 
78 genes from 17 loci that had been found by previous GWAS meta-analyses to have a significant 
association to one of the six traits. In addition, we had access to the first five principal compo¬ 
nents of genome-wide genotypes. The goal in Service et al. ( 2014[ > is to identify which variants in 
these loci are most likely to directly influence the observed variability in the three distinct lipid 
traits. 

Data cleansing and filtering are described in detail in the supplement; here we limit our¬ 
selves to note that for the purpose of the simulation study, the collection of variants was pruned 
to eliminate 550 variants observed only once and to obtain a set of variants with maximal pair¬ 
wise correlation equal to 0.3 by removing another 558 variants. We excluded singletons from 
consideration since it would not be possible to make inference on their effect without strong as¬ 
sumptions. Multiple considerations motivated our choice of selecting a subset with only modest 
correlations: (a) correlated variants make the convergence of MGMG problematic, which might 
impair our ability to understand the inference derived from the posterior distribution; more 
importantly, (b) it is very difficult to evaluate and compare the performance of model selection 
methods in the presence of a high correlation between variants; and finally, (c) statistical methods 
cannot really choose between highly correlated variants and the selection among these needs to 
rely on experimental studies. Let us expand on these last two points. Procedures that build a 
multivariate linear model, such as the Lasso, would select one out of multiple highly correlated 
variants that have some explanatory power for the response; approaches such as BH marginal 
would instead tend to select them all; and Bayesian posterior probabilities for each of the variants 
would reflect the fact that substitutes are available: there will be multiple variants with elevated 
(if not very high in absolute terms) posterior probability. 

It becomes difficult to meaningfully compare FDR and power across these methods, sub¬ 
stantially reflecting the fact that the problem is somewhat ill-posed: if multiple highly correlated 
variants are available, any of them can stand for the others, and it is arbitrary to decide on purely 
statistical grounds that one belongs to the model while the others do not. Since our goal here is 
to understand the operating characteristics of the procedures, we found it useful to analyze them 
in a context where the target is well-identified and the results easily comparable. 

After the described pruning, the genetic data used in the simulations contains 5335 subjects 
and 768 variants. Genotypes were coded with minor allele counts, and missing values (0.04% of 
genotypes) were imputed using variant average counts for consistency with previous analysis. 
Observed minor allele frequencies range from 2 x 10^^ to 0.5, with a median of 0.0009 and a 
mean of 0.02. There are 628 variants with MAF<0.01. Annotation information was obtained as 
in [Service et al. ( |2014| >, resulting in 61% coding, 34% UTR, and the remainder intragenic. Prior to 
applying the selection methods, the five genetic principal components along with the intercept 
were regressed out of both X and Y, and the columns of both were then standardized. 

When studying a real dataset, however, investigators might not be comfortable with such 
a stringent level of pruning; one might be concerned that variants with important effect are 
eliminated and that one is essentially reducing the information content of the sample. Indeed, 
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when analyzing real data, we used a much more comprehensive approach, as described in the 
case study section. 

3.4 Simulation scenarios 

We constructed two simulation scenarios: one to simply illustrate the advantages of the pro¬ 
posed priors and the other to investigate their potential in a set-up that models a real genetic 
investigation. 

Illustrative example: orthogonal X. 

We set n = 5000, p = 50, q = 5, and X = Ip ■■ ■ Ipf so that X^X = (n - l)7p. In 

generating and the responses, we want to cover a range of different signal-to-noise ratios. To 
achieve this, we sample values of the parameters using the distributional assumptions that we 
described in the specification of the priors. To explore the performance of the Across Traits and 
Across Sites models—both when they provide an accurate description of reality as well as when 
they do not—we use three rules to generate the probability with which each variant is associated 
to each trait, (a) We sample one sparsity parameter oo for each trait and keep it constant across 
variants, (b) We sample a probability for each variant and keep it constant across traits. Finally, 
(c) we define groups of five variants and sample one probability Vg of causality for each group 
of variants and each trait. Rules (a)-(c) are most closely reflected in the prior structure of the 
basic. Across Traits and Across Sites models, respectively; and we indicate them as exchangeable 
variants, pleiotropy, and gene effect. We generate 100 datasets per rule, each with q responses, 
and analyze them with the described set of approaches. When using Bayesian methods, we rely 
on non-informative priors (see supplement for details). 

Actual genotypes X. 

To explore the potential power and FDR in the analysis of the dataset with three lipid traits, we 
generate artificial phenotypes starting from the available pruned genotypes. We consider a mix¬ 
ture of possible genetic architectures. In the construction of each dataset, (a) one gene is selected 
uniformly at random for each phenotype and 3^ of its rare variants are causal (gene effect); (b) 
40 distinct common variants are selected uniformly at random and each has probability equal 
to 0.1 to be causal for each of the phenotypes (thereby substantially representing trait-specific 
variants); and, finally, (c) 10 additional common variants are selected uniformly at random and 
each has a probability 0.9 to be causal for each phenotype (pleiotropic effects). This results in 
traits that are on average determined by 3-4 rare variants in one gene, 4 common variants with 
effects on one trait only, and 9 common variants with effects across multiple traits. We generated 
a total of 100 such datasets, as detailed in the supplementary material. 

3.5 Data Availability 

The sequencing and phenotype data are available on dbGaP. The Northern Finland Birth Co¬ 
hort 1966 (NFBC1966) study accession number is phs000276.v2.pl. The Finland-United States 
Investigation of NIDDM Genetics (FUSION) study accession number is phs000867.vl.pl, with 
the sequencing data in substudy phs000702.vl.pl. In both cases, the sequencing data used in this 
paper has molecular data type equal to Targeted Genome' rather than 'Whole Exome'. 
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4 Results 


4.1 Simulations 

Illustrative example. 


BH marginal 

Lasso - 

Across Traits —•— 

Unadjusted 

BHfull —^ 

— basic —•— 

Across Sites 




FDR 

Figure 4: Power (top) and BFDR (bottom) as a function of empirical FDR in the illustrative 
example. Each color indicates a different variant selection approach (see legend at the top). 
Displays in different columns are for different data-generating mechanisms. 

Figure showcases the possible advantages of the priors we have described. The plots on 
the top row compare the empirical FDR and power of the different variant selection methods 
on the datasets with orthogonal X. Points along the curves are obtained by varying tuning 
parameters and averaging the resulting PDF and power across 100 simulated datasets. Our 
setting is such that BH full, BH marginal. Lasso and the basic Bayes model have very similar 
behaviors: the Across Traits and Unadjusted models achieve the highest power per FDR in the 
presence of pleiotropy and the worst power per FDR in the presence of gene effects; in contrast, 
the Across Sites model has maximal power in the presence of gene effects and worse power in 
the presence of pleiotropy. While it is not surprising that the most effective prior is the one that 
matches more closely the structural characteristics of the data, it is of note that the loss of power 
deriving from an incorrect choice of the Across Traits or the Across Sites model is minimal for 
FDR values lower than 0.2, which are arguably the range scientists might consider acceptable 
(see supplementary information for a detail of these values). On the bottom row of Figure]^ we 
compare the estimated BFDR with actual FDR for the Bayesian models; here the most serious 
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mistake is in underestimating FDR, which would lead to an anti-conservative model selection. 
Once again it can be seen that the best performance is obtained with the prior that matches 
the data-generatmg process. Besides this, it is useful to analyze the behavior of the Unadjusted 
approach: its power increase per FDR in the presence of pleiotropy is less pronounced than that 
of the Across Traits model, substantially because the Unadjusted approach is too liberal, with a 
BFDR which is significantly underestimated. This is in agreement with the lack of adjustment 
for multiplicity indicated by (|^. Results for alternate hyperparameters are in the supplement. 

Generating phenotypes from actual genotype data. 

Figure shows the performance of the variant selection methods in the analysis of traits gen¬ 
erated from actual genotype data, further emphasizing the potential gains associated with the 
proposed strategies. For given FDR, both the Across Traits and Across Sites priors lead to an in- 

BH marginal Lasso - Across Traits — •— 

BH full —•— basic — •— Across Sites 




FDR 

Figure 5: In the top portion, power and BFDR as a function of empirical FDR in the simulation 
from actual genotype data. In the lower panels, power is calculated separately for the different 
types of variants: rare, common with trait-specific effects, and common with pleiotropic effects. 
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Table 1: FDR and power for specific choices of selection parameters applied to simulated traits 


Variant selection criteria 

FDR 

power 

BH full with a = 0.05 

0.037 

0.32 

BH marginal with a = 0.05 

0.079 

0.37 

Lasso min error A 

0.810 

0.63 

Lasso 1-se A 

0.046 

0.26 

basic BFDR < 0.05 

0.039 

0.37 

Across Traits BFDR < 0.05 

0.054 

0.45 

Across Sites BFDR < 0.05 

0.051 

0.40 


crease in power over the other methods. This is due to the fact that phenotypes are generated 
assuming both pleiotropy and contributions from multiple rare variants in the same gene (gene 
effects). In the lower portion of Figure]^ we separate the power to recover rare variants with gene 
effects from that for trait-specific common variants (a; = 0.1) and from that for common variants 
with pleiotropic effects (a; = 0.9). As expected, the gains of Across Traits and Across Sites are for 
the portion of genetic architecture that is accurately reflected in these priors. The estimates BFDR 
are accurate, indicating that all three Bayesian priors correctly learned r and the probabilities of 
function. 

Finally, while we have relied on ROC-like curves to compare different approaches as the 
value of their tuning parameters vary, it is useful to focus on the operating characteristics of the 
standard ways of selecting the tuning parameters. By convention, the target FDR for BH is usu¬ 
ally 0.05. For Lasso selection, the function cv. glmnet provides two choices for A: minimizing 
the cross-validation error and using the one-standard error rule. In Bayesian approaches, one 
can select variants so that global BFDR < 0.05. Table [l compares FDR and power for these selec¬ 
tion parameters; the Bayesian methods appear to control the target FDR and arguably result in 
better power. Analogous summaries for other decision rules are included in the supplementary 
material; here we simply remark that including variants such that BFDR < 0.05 in this dataset 
was practically equivalent to selecting variants with posterior probability larger than 0.7. We will 
capitalize on this observation for the real data analysis. 

The supplement details the results of another set of simulations along the lines of a tradi¬ 
tional investigation of pleiotropy versus coincident linkage; we give a very brief summary here. 
In the case of separate causal variants, the Across Traits prior may have a slight loss of power but 
is still much better than BH with p-values from the full model. In the case of pleiotropy, however, 
the Across Traits prior clearly has greater power per FDR. 


4.2 Case Study: the influence of 17 genomic loci on lipid traits 

We now turn to the analysis of the three lipid traits in the Finnish dataset. While resequencing 
data comes from 17 loci identified via GWAS, prior evidence of association is available only be¬ 
tween some of these loci and some traits. In particular, four loci have no documented association 
with any of the three lipid traits we study; we include variants from these loci in the analysis as 
negative controls. (This is different from the work in Service et al. ( |2014 |, which examines only 
variants in loci specifically associated with each trait.) 
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Table 2: Summary of selections for BH with a = 0.05, Lasso with A chosen by cv.glmnet, 
and Bayesian approaches with ^ = 0.7. The columns labeled R and V* give, respectively, the 
number of variants selected across the entire study and in the four loci with no prior evidence 
of association to any of the lipid traits analyzed (CRY2, G6PC2, MTNRIB, and PANKl). BFDR 
reports, for Bayesian methods, the Bayesian FDR computed separately for each trait. 



HDL 

EDE 

TG 

Variant Selection 

R 

V* 

BFDR 

R 

V* 

BFDR 

R 

v* 

BFDR 

BH full p-values 

13 

0 


3 

0 


5 

0 


BH marginal p-values 

22 

0 


6 

0 


10 

0 


Easso min error A 

134 

12 


40 

4 


80 

6 


Lasso 1-se A 

16 

0 


0 

0 


4 

0 


basic 

21 

0 

0.046 

4 

0 

0.012 

8 

0 

0.011 

Across Traits 

19 

0 

0.084 

8 

0 

0.059 

9 

0 

0.029 

Across Sites 

25 

0 

0.064 

5 

0 

0.063 

8 

0 

0.007 


Service et fl/.|(201^ relied on univariate regression to test the association between each trait 


and each variant with MAF>0.01 and on burden tests to evaluate the role of nonsynonymous 
rare variants. Bogdan et al. ( |rn press I re-analyzed the data relative to HDL with a set of model 
selection approaches, including the novel methodology SLOPE; to facilitate comparison with 
their results, we add SLOPE to the analysis methods considered so far. Groups for the Across 
Sites model were defined so as to mimic the burden tests in Service et al. ( 2014| >, which means a 
group with more than one variant contains all nonsynonymous variants with MAE<0.01 in the 
same gene. 

We start by analyzing the pruned subset of variants used in the simulation studies and 
postpone a more exhaustive search, noting again that this allows for a more straightforward 
comparison of the variants selected by different methodologies. Table compares the number 
of variants selected by various methods with specified tuning parameters. The column labeled 
V* shows the number of selected variants that are in a locus lacking any prior evidence of 
association to lipid traits. The Easso with A chosen to minimize cross-validated prediction error 
clearly results in far too many selections, so we discard this approach for the remaining results. 
Eor Bayesian approaches, the threshold ^ = 0.7 results in average BEDR approximately controlled 
at the 0.05 level. 

Eigure illustrates the model selection results for HDE. (Analogous displays for the other 
two phenotypes are available in the supplementary material, as well as a table detailing differ¬ 
ences in selections between approaches.) Each display corresponds to a locus, with turquoise 
shading (rather than orange) used to indicate prior evidence of association to HDE. Variants are 
arranged according to their genomic positions in the loci, and the values of their estimated co¬ 
efficients are plotted on the y-axis; with the exception of marginal BH, we display only nonzero 
coefficients. When available, a vertical black line indicates the position of the SNPs originally 
used to select the locus ('Array SNP'). 

There is substantial overlap among the results of various methods. Model selection ap¬ 
proaches seem to generally agree with the findings in Service et al. (2014 1 (with Easso 1-se the 
most discrepant, missing a number of the associations identified in Service et al. ( |2014 1; see sup¬ 
plementary material). Still, we can point to some significant differences. With the Across Traits 
approach we select two variants in two loci where no other method identifies any signal: in 
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Figure 6: Estimated variant effects on HDL. Each panel corresponds to a locus, the z-axis 
indicates the variant's genomic position and the y-axis its regression coefficient (with the excep¬ 
tion of BH marginal, only nonzero coefficients are represented). The color code of the panel 
titles indicates the presence/absence of prior evidence of association between the locus and HDL 
(turquoise/orange, respectively). Model selection methods are distinguished using plotting sym¬ 
bols, as indicated in the legend at the top. 
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CELSR2 and FADSl. These two loci have prior evidence of association to LDL and to all three 
lipid traits, respectively, and the Across Traits approach identifies pleiotropic effects. In contrast, 
the Across Traits approach does not select four very rare (MAF<0.001) variants considered rele¬ 
vant by more than one alternative method. While we do not know where the truth lies at this 
point, it is very hard to evaluate the effect of a rare variant on purely statistical grounds, and the 
outcome of the Across Traits model might well be the more reliable. 

The Across Sites approach identifies four variants that other approaches overlook. Three are 
rare variants in ABCAl: two missense rare (0.01>MAF>0.001) and one nonsense very rare (MAF 
0.00016); their discovery is facilitated by the fact that they are included in a group with multiple 
other significant variants. The fourth is a common variant in the MVK locus, for which there 
is prior evidence of association to HDF. Other approaches do not recover this simply because 
the signal in the locus is, in the dataset analyzed here, barely below the detection threshold; 
Across Sites has a slight advantage over the other Bayesian methods because grouping reduces 
the number of comparisons to account for. We note that ABCAl is a gene in which rare variants 
were found to have a role by the burden test analysis in Service et al. ( |2014| |. 

We have relied on the pruned dataset since low correlation across variants greatly facilitates 
the comparison of the selection results by different methods. However, mindful of the concerns 
of scientists unwilling to pre-screen the genotypic information, we have also carried out a more 
comprehensive analysis of this dataset, showcasing that this is indeed an option available to 
researchers. Details for these analyses are in the supplement, but we summarize them here. 
First, we have compared the results of four different levels of pruning (correlation less than 0.3, 
0.5, 0.7 and 0.9). We have found that very few values of Z^t change by more than 0.05 when 
different levels of pruning are used and that less stringent levels of pruning do not lead to 
substantially more findings—unlike when applying BH to marginal p-values. In fact, there is 
a greater tendency for variants to drop out of the selection set as correlated variants are added 
to X. Second, to completely eliminate pruning, we analyzed all variants in one locus along the 
lines of Servin and Stephens (2007[|; Hormozdiari et al. (|2014||; Kichaev et al. (|2014||; Chen et al. 


( |2015 l, using the basic prior and assuming that the number of significant regressors |Z| is no 
greater than five or six (depending on the total number of possible variants in the locus). We 
have restricted our attention to two loci only, those that showed stronger evidence of influencing 
HDF via multiple variants. The supplementary material offers a precise comparison of results, 
but it suffices here to note that the set Z of variants with the largest posterior density is equal to 
the variants selected among the original pruned set (correlation less than 0.3) by the basic prior 
for one locus and, for the other locus, the two sets substantially overlap. 


5 Discussion 

As the genetics community devotes increasing effort to follow up GWAS hits with resequencing, 
a number of suggestions on how to prioritize variants have emerged. In truth, while dealing 
with the same broad scientific goals, many of these contributions address different aspects of 
the problem and therefore should be seen as complementary rather than alternatives; taken 
together they provide scientists with useful suggestions. Annotation information has been shown 
to be quite useful when the set of variants under consideration is sufficiently diverse. It is 
important to account for the correlation across variants to avoid paying attention to SNPs that 
are only 'guilty by association.' Bayes's theorem and Bayesian statistics are a natural way of 
dealing with the decision of which variants to pursue. In this context, others have studied 
(a) how to choose priors that incorporate annotation information, tuning their parameters with 
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available datasets; (b) how to approximate posterior distributions of variant effects; (c) how to 
sample from the posterior distribution using efficient MCMC or variational schemes; and (d) how 
to efficiently evaluate posterior probabilities for a set of variants. Here we focus on another 
aspect of prior selection: describing how partial exchangeability assumptions can be used to 
borrow information across traits and neighboring sites, while maintaining an effective control for 
multiplicity across variants and fitting multivariate regression models that estimate the specific 
contribution of each associated site, while accounting for others. We briefly refer to some of the 
most direct antecedents of our proposed priors to underscore relevant differences. 

Yi et al. (2011| proposed the use of hierarchical priors to capture effects of rare variants 


through groups, similar to the Across Sites model. However, their proposal does not incorporate 
sparsity considerations, resulting in the estimate of a nonzero effect for each variant and each 
group and therefore not engaging in model selection. Quintana et al. ( |2011 2012[ ) took an addi¬ 
tional step towards the Across Sites model by incorporating sparsity via the indicator variable Z. 
They considered only rare variants, used the same effect size for all rare variants in a genomic 
region, used the MLE for the effect sizes rather than integrating them out, and, most importantly, 
controlled sparsity by using A^ = 1 and Eg = pg in the prior for Vg rather than introducing 
another layer of indicator variables in the hierarchical prior—all of which means their approach 
has less flexibility and less learning. 

The Across Sites prior also echoes the proposal of Zhou et al. ( 2010[ | who suggested the 
use of group penalization in Lasso to estimate multivariate sparse models while encouraging 
coordinated selection of rare variants in the same gene. This computationally appealing approach 
has not become as popular in genetics as in many other fields, possibly because of the difficulties 
connected with the selection of its tuning parameters when model selection is the goal. Cross 
validation is often used to determine the appropriate level of penalization; while this works well 
for prediction purposes, its performance is less than satisfactory in identifying variants with truly 
nonzero coefficients (as illustrated by our case study). [Alexander and Lange ([2011 >,[Valdar et oT. 


(20121, and, most recently, Sabourin et al. ( 2015| explore coupling resampling techniques with 
Lasso penalization to improve model selection. This not only increases computational costs but 
also greatly reduces the initial simplicity of the model. As documente d in [Bogdan et al. ( in press) , 
identifying a single A value that performs LDR control is challenging; Yi et al. ( [2015[ l irivestigates 
this task in the context of GWAS and provides guidelines. The final model selection of these 
machine learning approaches uses complex rules; in contrast, the Bayesian models we described 
are based on easy to interpret parameters. 

The use of hierarchical Bayesian methods has ample precedents in eQTL studies, where 
they have been used to correct for multiplicity ( Kendziorski et flr][2006[ l and to increase power of 
detecting variants affecting multiple traits. In our presentation of the Unadjusted and Across Traits 
approaches, we referred to methods proposed by[Jia and Xu[([2007[| and[Bottolo et fll.[([2011[|. More 
recent work ( [Llutre et a?. [2013 Li et fll.[[2013[ l has focused on the identification of local (cis) effects 
across tissue, and considered models with only one functional variant. The recent contribution 
by Chung et al. (20141—^which appeared while this work was in preparation—underscores as we 
do the importance of learning both across sites and across traits to prioritize variants. These 
authors, however, work with p-values from GWAS studies, rather than actual resequencmg data. 

Having clarified the scope of our contribution, we want to briefly mention how it could be 
extended and combined with suggestions by others. Lirst, let us point out that while in the 
simulations and in the analytical approximations we assumed n > p, this restriction is by no 
means necessary to the Bayesian model we describe. On the contrary, the priors we propose— 
by learning sparsity and giving positive probabilities to configurations with some j6t, = 0—are 
well-suited to the case n < p. The real challenge in dealing with GWAS-type data would be from 
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a computational standpoint: increased mixing for MCMC as described in Xu et al] (|2014| and 


Guan and Stephens ( 2011| | or other algorithmic improvements ( [Carbonetto and Stephens 
Hormozdiari et fl?.||2014[| would make our approach more widely applicable. 


2012 


Another extension that is easily achieved is the combination of the Across Traits and Across 
Sites priors. Most immediately, the group indicators G in Figure]^ can be made trait-specific and 
linked across phenotypes with the same approach used to link the in Figure 

It is certainly possible to combine the partial exchangeability aspects of our models with a 
prior that incorporates annotation information. Refer, for example, to the Across Traits prior in 
Figure]^ Gurrently, the distribution on Wi,..., Wp, indicators of functionality of the variants, is a 
beta-binomial. However, it is trivial to change it to a mixture of independent logits, with the linear 
model component including an intercept effect—^which would capture the overall sparsity—and 


a linear combination of annotation indicators (Veyrieras et al. 
2014l l. 


2008 Kichaev et al. 2014 Pickrell 


Since our focus has been on specification of the prior, we have not paid much attention to the 
data-generating model, which could certainly be improved. Specifically, we want to underscore 
the fact that using a mixed-model approach might be advisable to account for population struc¬ 
ture (|Kang et flZ.|2010]l and when analyzing many phenotypes whose quantitative value might be 
influenced by confounders ( [Zhou and Stephens|2014 1 or simply by genetic variants not included 
in the model. 

In conclusion, we want to emphasize the increasing importance in human genetics of models 
that account for pleiotropy. 'Big data' in genetics has often been equated with the abundance of 
sequences, and these certainly pose a number of management and interpretation challenges. 
Our increased acquisition capacity will also result, however, in the collection of a large number 
of phenotypes; gene expression, MRI scans, and mass spectrometry are just some examples of 
the large-scale phenotypmg efforts underway. Now that DNA variation has been extensively 
described, annotating this appears as a fundamental challenge; the rich phenotypic collections 
increasingly available have a major role to play. After all, what better way of establishing if 
a variant has some functional impact than looking for its association with any trait available? 
Bayesian models that allow one to estimate the probability with which a variant has functional 
effects across phenotypes are likely to be very useful. In this paper, we have described a first step 
in this direction. 
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A Mathematical details for the basic prior 


First we integrate j8 out of the posterior distribution as in ( |Chen et~ar. |2015| |. Since y = ^ 

with {^^\Z,p,T) ~ M (^0, yEz) and e ~ A/'(0, jin), 

{y\Z,p,T) (^0,h„ + jXz^zXiy ( 6 ) 


While the likelihood can be written directly from this, a few manipulations give it in a more 
convenient form. Sylvester's determinant theorem (|Harvill^|2008 p. 420) implies 


- t-2|Z| 


. 2 |z| det(Ez) 


det(J„ + t'XzEzX^) = det(J|z| + r^EzX^Xz) = r . 

and a generalization of the Woodbury matrix identity (|Harvill^|2008[ p. 428) implies 


{In + t'XzEzX^)-! = In- Xz{t-^'^z' + ^Z^z)~'X^z = hr - XzOzX 


2v-l 




z- 


Consequently, Pr(y|Z,p, t) = rl^rdeFE For the null model Z = 0, ^ shows 

that the covariance matrix is p^^In, so in this case S| = y^y and the ratio of determinants is set 
equal to one. 

Next multiply Pr(y|Z,p, t) by the prior density function for p and integrate to obtain 


det(Qz)^^^ 

det(Ez)^'^^ 





det(Qz)^^^ 

det(Ez)^/^ 




since the integrand is the density function of Gamma (ctp + |, Ap + jSl), up to a normalizing 
factor. Hence, the marginal posterior for Z and t is given by (|^. 

Along the lines of Malsiner-Walli and Wagner (2011|l, we present an approximation of the 
posterior expected value E[Zi,|Z[_j,],T,y]. If Z and Z are equal except that = 0 and Z^, = 1 for 
one V, then 


fz,r{tT\y) ^ Fr(Z^ = l|Z[_^],T,y) ^ E[Z^|Z[_^],T,y] 

fz,r{Z,T\y) Pr(Z^ = 0|Z[^„],T,y) 1 - E[Zj,|Z[_^],T,y]' 


We will use several assumptions in order to simplify the expression on the left and then solve 
for E[Zi,|Z[_j,], T,y]. Consider the case when the columns of X are orthogonal, which implies 
X^X « nip and the two choices of Ez are essentially the same. Furthermore, {x^,, y) = (x^,, X/S + 
e) ~ njSv + {xn,e), which is distributed as Af{0, ^{Z^nr^ + 1)); in this context, distinguishing 

signal from noise requires nr^ » 1, so we assume this to be the case. Consequently, ~ 
(n + T^^)Ip « nip, which in turn implies S| is approximately equal to the residual sum of 
squares (RSS) for the model indicated by Z. Finally, reflecting the results of current GWAS, we 
assume that the portion of variance explained (PVE) by the loci in consideration is rather small, 
so RSS is not much less than y^y « n for any model. If one further chooses iXp = \p j, then 


fz,r{^,T\y) ^ 

1 /z(2) 

/det(Q2)det(Ez)y/' , 


fz,r{^,T\y) ' 

~ T/z(Z) 

ivdet(E 2 ) det(Qz)y ' 
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Properties of the beta and gamma functions give that the ratio of fz values is + |Z|)/ (B^ + 
p — |Z| — 1). Furthermore, Y/det(Q 2 ) det(Ez)/ deiiJL^) det(Qz) ~ 1/ 'Jn. Finally, S| ~ y^y — 
^E»:Z„=l(4y)^ so 

S|- . (xjy)2 _ 2 

S| ^ S| ^ ny^y “ 

Substituting these results into 0 and (|^ gives (0. 


B Mathematical details for learning across traits 


While the Unadjusted prior is not useful, we include its marginal posterior density here for com¬ 
pleteness. Its derivation is very similar to that of the basic model, so we focus on the differ¬ 
ences. A priori the rows of Z are independent and each has a beta-binomial distribution, so 
/z(Z) = Ylv=i B(Ai, + Sv,Bi, + q — Sv)/B{Az,,Bz,). Furthermore, the columns of Y are independent 
given Z, )3, and p, and similarly for the columns of j8; so 


/z,t(Z,t|Y) oc/.r(T)/z(Z) 


n 


f=l 



( 2 +“p ) 


det(QzJ^^^ 
rlZfl det(Ez,)^/^ 


(9) 


If Z and Z are equal except that Z^,t = 0 and Z^,t = 1 for one v and one t, then the same 
calculations as for the basic model give that (0 simplifies as 

/z,t(Z, T|y) ^ ^ _ A^, + Sy , ^ 2 \-n/2 

fz,T{^,T\y) ^ T^/nBj, + q-s^,-l^ 

This leads to the approximation (0. 

Next we consider the Across Traits prior. The posterior density is the same as in (0 except 
that fz is replaced by 

/w,z(W,Z) = J Fr{Z\W,v)Fr{W\cow)Pr{v)Fr{cow)dvdcow 

_ B(Avv + I W|, Bw -|- p — I W|) .p-j- B(Ai,-|- Sy, Bj, -|-— Sd) 

B(Aw, Bw) i;:W„ = l B(Ai„Bi,) 

provided that Zj^t = 0 for all v such that W^, = 0—otherwise, /w,z(W,Z) = 0. 

To derive the approximation for the posterior expected values, consider W and W that are 
equal except that FIv = 0 and = 1 for one v. Choose Z consistent with W, which means 
/w,z(W,Z) / 0. Since we are using a subscript to denote a column of a matrix, we use a 
superscript as in Z^ to denote row v of Z. Furthermore, zt“^l denotes the sub-matrix of Z 
obtained by deleting row v. Straightforward modification of 0 gives 

Ez/w,z,.(W,Z,t|Y) _ E[W,|W[„,],zM,t,Y] 

/w,z,t(W,Z,t|Y) 1-E[W,|W[_,],ZM,t,Y]' 

where the summation is over all Z such that Zt“^l = Z[“^l. Furthermore, for any such Z, 

/w,z,t(W, Z, t| Y) ^ A I Wj_j,j I BjAjj-\-Sjj,Bi,-\-q s^) .p-i- 1 , 2\—n/2 

/w,z,t(W,Z,t|Y) ^ Bw + p-|W[_„]|-l B(A„,B,) ^ ~' 
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Hence, E[W^|W[_^],Z[-^U,Y]-i « 1 + Y) with 


i2,(W[_„],T,Y) 


A]/m -\- 



B(Ai, -|- 

Z^ 

,Bz)-\-c] — 


Bw + P — 


_ 1 Zj 

^ Z'' 

B(A 


1 1 i Vvt) ' 


where the summation is over all 2*? possible values of Z^. 


C Mathematical details for learning across sites 

The derivation of the marginal posterior for learning Across Sites consists of straightforward 
modifications of previous calculations. The marginal posterior distribution is as in except that 
the prior /z(Z) is replaced by the joint prior for Z and G, which is 

f (r 7 ) + |G|/Sg + - |G|) p|- B(Ag + Sg,Bg + Pg-Sg) 

-- B(A^) -Jj, — mAy — 

pg>^ 

provided that Zj, = 0 if = 0, and Zj, = 1 if = 1 and = 1—otherwise, fG,z {^'= 
0. Similar to the preceding approximations, E[Gg|G[„^],Z[_^],T, y] 1 l + l// 2 g(G[^g],T,y), 
where for Pg > 1 


l2g(G[_g],T,y) 


Z.G + 

G[-dl V- 

B(Ag + 

^g 

'^g~^Pg 

Z,l) 

Bc + r- 

^[-g] 

-l^—i 

B(Ag,Bg)(Tv/n)l^^l 


^■7(A=S, 

Zv=l 


with the summation being over all possible values of Zg. When pg = 1, however, the sum¬ 
mation is replaced by (1 — (T-\/n), where v is the lone variant in group g. Hence, the 

posterior conditional probability of Gg depends upon the overall number of groups and the num¬ 
ber of groups considered relevant, while the posterior conditional probability of Z^, given that 
Gg = 1 depends on the number Sg of variants in the same group that are deemed functional. 
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Supplementary Material 


Choice of hyperparameters 

We state here the hyperparameters used in most of our computations; the effect of alternate 
hyperparameters is discussed below. For fine-mapping, there is considerable uncertainty about 
the likely value of cu, so we use the uninformative prior with A^j = = 1. Since n > 5000 in 

our examples, we use = Ap = 10 to ensure that a.p = \p ^ 

The bounds on the uniform prior for t justify additional explanation. The lower bound 
could be zero, but this may lead to increased computational times because very small t implies 
small nonzero values in jS are more likely, which may make the MCMC iteration accept larger 
models. The assumption 1 suggests choosing Ti ~ 1/ \^n instead. Furthermore, T^/p is 

approximately the variance of the nonzero coefficients since the diagonal entries of either Zz are 
approximately one. The standardization of y and the assumption of small PVE implies p ~ 1. 
Consequently, T 2 should be at least three times as large as the largest likely value of \f5z,\. Either 
marginal or multivariate linear regression of our actual data gives all coefficients with absolute 
value less than 0.2, so we use Ti = 0.01 and T 2 = 10 in our analyses. 

MCMC sampling 

We use Gibbs sampling to alternate between updating t and updating the indicator vari¬ 
ables. We use Metropolis-Hastings for each update. We only update the indicator variables 
once for each update of t, but it might be more efficient to try multiple indicator variable jumps 
between each jump of r. 

To define the notation, the Metropolis-Hastings algorithm to generate samples with density 
proportional to /(•) draws a proposal u from a jumping distribution Q(m|m(’”)), where is the 
previous sample. Next, compute 


/(fl) 

Q /(«('")) ■ 


( 10 ) 


The new sample is u with probability min(R, 1) and is otherwise. To avoid overflow, we 
actually compute logR and then take its exponential. The remainder of this section gives the 
definition of the jumping distributions used. 

The jumping distribution for f given t is Af{T, (t 2 — Ti)^/ 16 ) but truncated to the interval 
(h/Tz)- When |Z| = 0 , the marginal posterior density is independent of t in (ti, T 2). Hence, if 
all effect sizes are very small, the t samples may be essentially uniform. This is because the data 
do not enable learning the effect sizes and has nothing to do with the MCMC algorithm or its 
convergence. 

With one level of indicator variables, our proposal distribution Q(Z|Z) is nonzero only if 
Z — Z = Se^, for <5 G {—1, -|-1} and one v, where is the vector with entry v equal to one and all 
others equal to zero. Furthermore, for such Z and Z, 


Qz(z|z) 


p 

P-|Z| 

l-^+ 

|z| 


if |Z| = 0 or |Z 
if ^ = +1 
if ^ = -1 


P 


( 11 ) 


where 7r+ is the probability of adding a variable to the model. 

Next we define Qwz(W/Z|W, Z) for learning Across Traits. The first step at each iteration 
is to choose S from {—1,0,-|-1} with probabilities Til, tTq, and tt^, respectively, with obvious 
modifications if, say, 3 = —1 is not possible because |W| = 0. Then W — W = 3ev, where 
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variant v is chosen uniformly from the possible candidates. The proposal Z will be the same as 
Z except for changes to Z^. If = —1, Z^ = 0. If ^ = 0, pT] | gives Q'^{Z'"\Z^) with q replacing 
p. If ^ = +1, then draw Z® based on the prior for (Z|W): first draw from the beta-binomial 
distribution 


p / + q Sp) f (j\ 

^ B{Ay,By) Uy 


( 12 ) 


and then sample Sy distinct entries uniformly from {1,... ,q}. The overall jumping probability is 
Qw,z(W,Z|W,Z) = Qo(Z|W,W,Z)Qw(W|W), where 


Qo(Z|W,W,Z) 


' 1 

< Q|(zlz") 

B(ApH~5 V Bp 


if = -1 
if = 0 

ifS = +l 


and the nonzero values of the mass function for (W|W) are 


Qw(W|W) = < 


Aw 

1 

P 


{l-nX)p 

K 

p-|W| 


|W| 


if^ = 0 
if |W| = 0 

if |W| = p and 3 = —1 
if 3 = -|-1 and |W| ^ 0 
if ^ = —1 and |W| ^ p 


(13) 


(14) 


The value of fiw is irrelevant because it will be the same for Qiv(W|W) and Qw(W|W), which 
will consequently cancel each other in p0| |. Furthermore, if 3 = +f, then Qo(Z|W, W,Z) in the 
denominator in p 0 | will cancel the leftover factor in fw,z{^>^) hr the numerator; whereas for 
^ = — 1 a factor in/w,z(W, Z) cancels Qo(Z|W, W,Z). 

The mass functions for the Across Sites prior are more complicated, but the fortunate cancel¬ 
lations still occur. First, pg replaces p in ([TT|) to give Q’^{Zg\Zg) when 3 = 0, and ([T^ becomes 


Qo(Z|G,G,Z) 


QI(Z,|Z,) 

1 

Pg^^g) 

b(a,,b,) 


if <5 = -1 
if ^ = 0 

if ^ = -|-1 and Pg = f 
if ^ = -|-1 and Pg > f 


(15) 


Let G* = 1 if = 1 and pg > 1; otherwise it is zero. Assuming that at least one group has 
Pg > 1, the nonzero values of the mass function for (G|G) are 


Ac 
1 
r 

n*_ 

(l-TT* )r 

r-|G| 

n*_ 


|G| 

(l-7r2)(r-|G|) 

n*_ 

(l-^o*)|G| 


if^i = 0 
if |G| = 0 

if |G| = r and 3 = —1 
if ^ = -t-1 and |G*| 7 ^ 0 
if 3 = -1,|G*| 0, and |G| / r 

if^ = +l,|G*| =0and |G| 7^0 
if 3 = -1,|G*| = 0, and |G| + r 
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Qg(G|G)= 


(16) 















In our experiments, 7r+ = 0.5 for the Z proposal, while tTq = 0.5 and n*_ = = 0.25 for 

the W or G proposal. One chain starts with no variants in the model while another starts with 
all variants in the model. For the other two chains, the variants are ranked by the absolute value 
of their correlation with yt, and the }t most strongly correlated have Z^t set to one. For the actual 
phenotypic data, }t is 10 for one chain and 20 for the other. For the simulations, one chain has 
}t equal to the expected number of variants that are causal for trait t (as determined by the co or 
V used to generate Zf), and the other chain has Jt double that. The burn-in interval was 10,000 
iterations. For the simulated data, 500,000 samples were then used to compute the averages, 
which took a couple of hours or less (typically much less, especially when using only one trait 
at a time) running the chains sequentially. For the actual data, 10 times as many samples were 
used for the averages; the chains were run in parallel and each MCMC sampling took 2-3 hr. We 
always used the g-prior. 

To assess convergence of the averages Z^t, we computed the average for each chain and then 
set AZjjt equal to the difference between the maximum and minimum averages over the four 
chains. For each prior applied to each set of 100 datasets, at least 95% of AZ^t were less than 0.05. 
For the priors that incorporate W or G, at least 95% of the analogous values of AW^, or AGg were 
less than 0.05 except that only 92% of AG^ satisfied this condition when the Across Sites prior was 
applied to the traits simulated from the genotype data. 

Evaluation of variable selection performance 

The p-values for use in the BH procedure are from ordinary least squares, using either one 
variant at a time ('marginal') or all variants simultaneously ('full'). For each dataset, we apply BH 
to all traits simultaneously, which means each BH test has pt] hypotheses. For the performance 
curves, we choose a target FDR, say a, and apply BH with that a to each dataset, computing 
the false discovery proportion (FDP) and of discoveries of true null hypotheses for each dataset. 
We then average over the 100 datasets to obtain the FDR and power for one point in our plots. 
For our performance curves, we used the FDR control targets {0.001,0.005,0.01,0.05,0.1,0.3,0.5}. 
Applying BH to the traits separately gave essentially the same results for the simulated data; the 
effect of testing traits separately in the actual data is discussed below. 

The performance curves for the Lasso show the effect of varying the penalty parameter A. 
We use the R function glmnet with its default parameters (with family=' gaussian' ) except 
without an intercept. This returns the fits for a set of A values, but that set depends upon the 
input data. We compute the FDP and power at each value of A output for that trait. Using 
locpolyO in the R package KernSmooth, we perform local linear regression on all these 
points for all datasets to approximate FDP and power as functions of A. When instead choosing 
a specific A by cross-validation, the one-standard error rule gives the largest value of A such that 
the cross-validation error is within one standard error of the minimum error. In this case, the 
choice depends upon the randomly drawn cross-validation bins, so for reproducibility we set the 
random seed to 1234 immediately before calling cv. glmnet. 

To compute the performance curves for Bayesian selection of variants, we choose threshold 
select using the rule Zj^t > ^ for each dataset, compute the FDP and power for each dataset 
(same as for BH), and then average over all 100 datasets to obtain one point in our plots. The 
values of ^ were 0.01£, 1 — 0.01£, and 0.1£ for £ = 1,...,9. 

Genotype and phenotype data 

Variants in dbSNP version 137 are indicated with their rs names, while other variants are 
referred to as v_c9_107548661, say, specifying chromosome and position (from GRCh37). Starting 
with the genotype and phenotype data in Service et al. ( 2014[ >, we removed 786 individuals that 
were missing values for any of the principal components or the three traits HDL, LDL, and 
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log-transformed TG, leaving 5335 subjects. From the sequencing data that had passed quality 
control, we removed two variants that were missing data for over 20% of the subjects as well as 
variants whose minor allele occurred only once in the remaining subjects, leaving 1326 variants. 

Additional filtering involved: (1) annotation, whose possible values were nonsense, mis- 
sense, coding synonymous, UTR 5', UTR 3', intron, and unknown; (2) PolyPhen2 predictions, 
whose possible values were probably damaging, possibly damaging, benign, and unknown; (3) 
marginal regression p-values from Service et fl?.| ( 201^ ; and (4) minor allele frequency (MAP). If 
all subjects have exactly the same genotype (including missing values) for two or more variants, 
we chose one of the variants using the ad hoc rules in Figure thereby eliminating 24 more 
variants. This X was used in the exact computation of the posterior probabilities in individual 
loci discussed briefly at the end of the Results section of the main paper. 

For the MCMC method of estimating the posterior probabilities across all loci simultane¬ 
ously, we further pruned to obtain a set of variants with maximal correlation less than a specified 
bound Cmax/ for reasons discussed in the main paper. Figure gives our algorithm to achieve 
this. The set Uq contains variants mentioned in Service et al. ( 2014| | and [Bogdan et al. ( in press) . 
It was introduced so that we only dropped such variants when they were correlated with each 
other, which simplified comparison with those papers. Furthermore, we manually chose which 
of two variants to drop in some cases, as we now explain. Service et al. (20141 list 39 associations 
between an 'Array SNP' and a trait. Of those, 16 had p-value less than 0.001 for their data and 
involved one of the three traits HDL, LDL or TG. Our filtering process dropped one of these: 
rsl0096633, which was an Array SNP for both HDL and TG, because it had correlation 0.96 with 


rs328, which has a nonsense mutation. Among the new discoveries reported by Service et al. 


( 2014[ ), rs651821 and rs2266788 in gene APOA5 are both associated with one of our traits, but we 
had to drop one of them because their correlation is 0.95. We dropped the latter because it had a 
slightly stronger correlation with the Array SNP rsl2805061. 

For almost all of our results, Cmax = 0.3. In this case, after filtering, we put rsl2805061 and 
rs3135506 (which is a missense mutation predicted to be probably damaging) back into X, which 
then had 768 variants. After imputation, all pairwise correlations in X are less than 0.3 except 
for two that are slightly greater. There are 628 variants with MAF<0.01. To explore the effects of 
pruning the variants, we also consider X obtained by setting Cmax to 0.5, 0.7 and 0.9, resulting in 
968,1042 and 1124 variants, respectively. 

Using the phenotype data (via the p-values) when filtering the variants could affect the 
accuracy of methods to estimate or control FDR, but our previous experience has been that this 
does not occur in this situation. 

Because substituting means for missing genotype values is not accepted practice, we com¬ 
ment on our use of this approach. Since we discard the two variants with very many missing 
values, only 0.04% of the values in X (with Cmax = 0.3) are missing, and only six variants are 
missing data for more than 1% of subjects—one variant is missing data for 2.2% of subjects and 
the next worst is missing 1.3%. With Cmax = 0.9, only 0.08% of the values in X are missing. 

Each of the traits (or its logarithm) were regressed on age, age^, and indicator variables for 
sex, oral contraceptive use, pregnancy status, and cohort. The five genetic principal components 
along with the intercept were regressed out of both X and the residuals Y from this regression, 
and the columns of both were then standardized. 

Simulation scenarios 

For the datasets with orthogonal X, we set n = 5000, p = 50, q = 5, and X = Ip ■ • • Ip) 


77pP P 

so that X^X = {n — 1)Ip. The columns of X were not centered because that would destroy the 
orthogonality. In generating phenotypic data, p = 1, t ~ Unif(0.045,0.063), nonzero ~ 
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For each of the 23 sets of variants for which all subjects have the same genotype 

If they are on the same chromosome and have the same annotation and PolyPhen2 prediction 
Choose one arbitrarily 

Else if exactly one variant has annotation 'missense' (none had 'nonsense') 

Choose it 

Else if no variant has annotation 'missense' 

Choose the one with annotation 'coding synonymous' (which meant discarding UTR 3' variant or 
variant without annotation) 

Else (this case only occurred once) 

Choose the one with PolyPhen2 prediction of 'probably damaging' (which meant discarding the variant 
predicted to be benign) 


Figure 7: Pseudocode for eliminating duplicate variants. 


For mode = 1,2 

For c = 0.9,0.8,..., Cmax 

variants = var.filter(variants, mode, c) 

Function var.filterjl/, m, cq) 

C is the absolute value of the correlation matrix of V, computed with the R option use="pairwise . complete . obs" 
Consider the entries of C,y in decreasing order 
If Cij < cq, return V 
y* = {vi,vj} 

If (m == 1) 

If Vi and Vj are on different chromosomes, continue 

Add to V* all on same chromosome as f,- and such that C;^ or Cji is greater than cq 
Uq contains any entries in V* that we particularly want to keep (see text) 

Lfi = argmin^, gy, p-value—but only consider common variants 

U 2 contains variants in V* with the 'worst' annotation, where the ordering is nonsense, then missense with 
probably damaging, then missense with possibly damaging, then missense with any prediction 
If Uq is not empty 

U is the entry in Ug with the highest priority 
Else if Ui is empty 
U = U2 

Else if U 2 is empty 
U = Ui 

Else 

Lf = (Ji n Lf 2 

If U is empty 
U = Ui 

If U is empty, U = V* 

Stop in the following rules when one gives single variant u: 

1. ueU 

2. u is variant in U with greatest MAF 

3. Choose u from U arbitrarily 
Drop u from V and C 


Figure 8: Pseudocode for filtering correlated variants. 
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N{Q,t^), and the distributions of the probabilities of association were as follows: 

Exchangeable variants one co ~ Beta(12,48) was drawn independently for each trait. 

Pleiotropy one (Uw ~ Beta(16,55) was drawn for each dataset, was drawn for each variant, 
and then ~ Beta (48,12) was drawn for each nonzero Wu. 

Gene effect coc ~ Beta (16,55) was drawn independently for each trait, Gg was drawn for each 
group of five consecutive variants, and then Vg ~ Beta (48,12) was drawn for each nonzero 
G,. 

These choices ensured 3> 1 so that signal was greater than noise, E[Zi,t] 0.2, and |Zt| 
exceeded 15 only about 5-15% of the time. Very roughly 


(yt/Yf) - +2(Xz,^Zt'^t) + ~ ~ + 1), 


so the generated traits also had variance close to one—^between 0.94 and 1.17. There were 100 
datasets, each with q = 5 traits. When applying the Across Sites prior to any of these data sets, 
the groups are the same as when generating the dataset with gene effect: the first group consists 
of the first five variants, the second group consists of the next five variants, and so forth. 

When simulating phenotypes from the actual genotype data, all rare (MAF<0.01) variants 
in the same gene with missense or nonsense mutations are in the same group. Every other 
variant—common or rare—is in a group by itself. Then Z was generated randomly as follows: 

• for each trait t, draw one group from amongst those with at least five variants and set 
Z-ut = 1 for 3^ rare variants (number drawn uniformly) in that group 

• draw 10 common variants and, for every trait t, set Z-ot = 1 with probability 0.9 

• draw 40 common variants and, for every trait t, set Z-ot = 1 with probability 0.1 

The rest of the process to generate each trait was the same as in the other simulations. Again, 
there were 100 datasets, each with q = 3 traits. There are 20 genes with at least five nons}mony- 
mous rare variants; 16 of them have fewer than 10 variants, but the others have 14-16 variants. 

We also created another set of simulations along the lines of a traditional investigation of 
pleiotropy versus coincident linkage, which we will call the pleiotropy-linkage experiment. We 
began with X created by setting Cmax = 0.7 in Figure We then took pairwise correlations 
between all variants with MAF >0.1. From the 10 pairs with absolute correlation greater than 
0.6 and in the same locus, we selected the five pairs listed in Table For each selected pair, we 
simulated 40 datasets, each with X restricted to the one locus and with two traits. In half the 
datasets, one of the variants is causal for both traits (pleiotropy, with each variant causal equally 
often) while the other datasets have one variant causal for one phenotype and the other variant 
causal for the other phenotype (coincident linkage). No other variants were causal. The rest of 
the process to generate each trait was the same as in the other simulations. 

Results for illustrative example 

For our Bayesian variable selection methods and several non-Bayesian methods applied to 
the datasets with orthogonal X, Figure [^compares the power for FDR less than 0.2. In this range, 
the power of the Across Sites prior applied to the dataset with pleiotropy is only slightly less than 
that of the basic prior or the non-Bayesian methods. Similarly, the loss of power of the Across 
Traits prior applied to the dataset exhibiting gene effects is slight when FDR is less than 0.2. 
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Table 3: Pairs of variants used in the pleiotropy-lmkage experiment. The total number of variants 
in the datasets for that pair is p. 


variants (MAP) 

locus 

P 

correlation 

rslll42 

(0.29) 

rs464218 

(0.45) 

CELSR2 

128 

0.70 

rs567243 

(0.49) 

rs560887 

(0.31) 

G6PC2 

62 

-0.66 

rsl5285 

(0.24) 

rs316 

(0.12) 

LPL 

25 

0.63 

rsl2445698 

(0.19) 

rs5801 

(0.18) 

CETP 

88 

0.65 

rsl801706 

(0.21) 

rs5882 

(0.39) 

CETP 

88 

0.66 


BH marginal 

Lasso 

Across Traits —•— 

Unadjusted 

BHfull — 

basic ♦ 

Across Sites —•— 




FDR 


Figure 9: Power (top) and BFDR (bottom) as a function of empirical FDR in the illustrative 
example, limiting to FDR<0.2. Each color indicates a different variable selection approach (see 
legend at the top). The different columns correspond to different data generating mechanisms. 
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BH marginal Lasso- Across Traits —•— Unadjusted —*■ 

BH full —•— basic ♦ Across Sites —•— 




Figure 10: Power (top) and BFDR (bottom) as a function of empirical FDR in the illustrative ex¬ 
ample with prior t ~ Unif(2,5). Each color indicates a different variable selection approach (see 
legend at the top). The different columns correspond to different data generating mechanisms. 


We also investigated the effect of using other hyperparameters. Figure shows the results 
when the prior for t was Unif(2,5) but everything else was the same as for the results in Fig¬ 
ure 1^ As suggested by the approximations given in the main paper for the conditional posterior 
expected values of the indicator variables, increasing t causes Z^t to decrease; so fewer variables 
are selected for given decreasing both FDR and power. This also means that Z^t underesti¬ 
mates the probability that variant v is causal for trait t, and BFDR greatly overestimates FDR. 
On the other hand, using the prior t ~ Unif(0.045,0.063) as was used in data generation gave 
essentially the same results as in Figure [^because the mean and variance of the r samples were 
usually less than 0.1 when using the uninformative prior r ~ Unif(0.01,10). 

Finally, we tried using different hyperparameters for the priors on the probabilities of causal¬ 
ity but with everything else the same as for the results in Figure]^ For the dataset with exchange¬ 
able variants, we again used uninformative priors for cuw and coq but all other Beta priors had 
parameters A = 12 and B = 48. For the dataset with pleiotropy, we changed only the hyperpa¬ 
rameters for cuw and v for the Across Traits prior, which now have the same values as during data 
generation, and likewise for the dataset with gene effects and the Across Sites prior. Figure [TT] 
shows that this did not have much effect on the results. This suggests there is not much advan¬ 
tage to using tighter priors, even if it were possible to guess the hyperparameters so accurately. 
Results for generating phenotypes from actual genotype data 
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BH marginal Lasso- Across Traits — •— Unadjusted — *■ 

BH full —•— basic —•— Across Sites —•— 



FDR 


Figure 11: Power (top) and BFDR (bottom) as a function of empirical FDR in the illustrative 
example with tighter priors for the probability of causality in some cases. Each color indicates a 
different variable selection approach (see legend at the top). The different colurrms correspond 
to different data generating mechanisms. 
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Table 4: FDR and power for selected values of ^ when applying Bayesian methods to simulated 


phenotypes with actual genoty 

oe data. 

FDR power 

FDR 

power 

FDR 

power 


? = 0.5 


0.7 


0.8 

basic 

0.100 0.41 

0.039 

0.37 

0.021 

0.35 

Across Traits 

0.133 0.51 

0.057 

0.45 

0.031 

0.42 

Across Sites 

0.143 0.46 

0.054 

0.41 

0.029 

0.38 


When applying the Across Sites prior to data with actual genotype data, the groups were the 
same as described above for simulating phenotypes from this data. This mimicked the burden 
Service ef'^ ( |2014| . 


tests m 


To create the performance curves in the lower row of plots in Figure 5 in the main paper, 
the FDR for each ^ is equal to that computed for the curves for all variants, but the power is 
computed separately for each type of variant. Table shows the FDR and power when using 
specified values of ^ to select variants after applying the Bayesian priors to the same set of 
simulations. The results when ^ = 0.7 are essentially the same as the results in the paper from 
controlling BFDR < 0.05. 

Next we discuss the results of the pleiotropy-linkage experiment. Over all the datasets, at 
least 95% of were less than 0.05 after 500,000 MCMC iterations, and when using the Across 


Traits prior at least 90% of AW^ were. Figure 12 shows the FDR and power for the basic and 


Across Traits priors as well as BH applied to both sets of p-values. In the case of separate causal 
variants, the Across Traits prior may have a slight loss of power when FDR is less than 0.5 but 
is still much better than BH with p-values from the full model. In the case of pleiotropic traits, 
however, the Across Traits prior clearly has greater power per FDR than the other methods, which 
have roughly the same power regardless of whether the traits have the same or different causal 
variants. For the variants other than the pair listed in Table 8% of the time the difference in 
Zyt between the two priors is less than 0.05. When the traits are pleiotropic, the differences are 
similar for the non-causal variant; but for the causal variant, Z^t for the non-causal variant is 
sometimes greater than 0.7 with the Across Traits prior even when it is less than 0.2 with the basic 
prior, consequently increasing power considerably. When the traits are not pleiotropic. Across 
Traits may inflate Z^f for the wrong variant of the pair, thereby increasing FDR slightly for given 
but very rarely increases Zz,t by much for the causal variant; so power is not increased. 

Results for the case study 

Figures 13 and |14| illustrate the model selection results for LDL and TG, respectively, which 
are interpreted in the same way as the figure for HDL in the main paper. Tables |5]-[^ list the 
methods that select each variant for association with HDL, LDL, and TG, respectively, for variants 
selected by at least one method. An x in a column indicates that the method selects that variant 
for that trait. The column labeled 'position' gives the position (from GRGh37) of the variant on 
its chromosome. 

For LDL and TG, no method identifies any rare variant, and the Across Sites selections always 
agree with the majority of approaches. The Across Traits approach selects a few variants that are 
noteworthy. For LDL, it selects a variant in each of three loci—GETP, GALNT2, and NGAN— 
where no other method identifies any signal. In each case, the variant that the Across Traits 
approach selects has a strong association with either HDL or TG. Two of these three loci have 
previous evidence of association with LDL. For TG, the Across Traits approach selects a variant 
in FADSl, where no other method selects any. This is the same variant that caught our attention 
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BH marginal basic • 

BH full • Across Traits • 



Figure 12: Power (top) and BFDR (bottom) as a function of empirical FDR in the experiment to 
check whether Across Traits can distinguish between pleiotropy (right column) and traits with dif¬ 
ferent correlated causal variants (left column). Each color indicates a different variable selection 
approach (see legend at the top). 
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Figure 13: Estimated variant effects on LDL. Each panel corresponds to a locus, the v-axis 
indicate the variant's genomic position and the y-axis its regression coefficient (with the excep¬ 
tion of BH marginal, only nonzero coefficients are represented). The color code of the panel 
titles indicates the presence/absence of prior evidence of association between the locus and LDL 
(turquoise/orange, respectively). Model selection methods are distinguished using plotting sym¬ 
bols, as indicated in the legend at the top. 
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BH marginal Lasso ^ basic Across Sites □ 

BH full + SLOPE V Across Traits O Array SNP I 
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Figure 14: Estimated variant effects on TG. Each panel corresponds to a locus, the x-axis indi¬ 
cate the variant's genomic position and the y-axis its regression coefficient (with the exception 
of BH marginal, only nonzero coefficients are represented). The color code of the panel ti¬ 
tles indicates the presence/absence of prior evidence of association between the locus and TG 
(turquoise/orange, respectively). Model selection methods are distinguished using plotting sym¬ 
bols, as indicated in the legend at the top. 
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Table 5: Comparison of selection results for HDL. 







BH 

BH 




Across 

Across 



variant 

position 

locus 

MAE 

full 

marg 

Lasso 

SLOPE 

basic 

Trait 

Sites 

PLOS 

v_ 

_c9_107548661 

107548661 

ABCAl 

0.0002 

X 

X 

X 

X 

X 


X 

X 

v_ 

_c9_107555091 

107555091 

ABCAl 

0.0005 


X 


X 

X 


X 

X 

v_ 

_c9_107555452 

107555452 

ABCAl 

0.0002 

X 

X 

X 

X 

X 

X 

X 

X 


rs76881554 

107578620 

ABCAl 

0.0043 







X 



rs2066715 

107588033 

ABCAl 

0.0551 


X 

X 

X 

X 

X 

X 

X 


rs2066718 

107589255 

ABCAl 

0.0148 

X 

X 


X 

X 

X 

X 

X 


rsl45183203 

107646756 

ABCAl 

0.0029 







X 



rs2575875 

107662494 

ABCAl 

0.2999 


X 

X 

X 

X 

X 

X 

X 

v_ 

_c9_107665945 

107665945 

ABCAl 

0.0002 







X 



rsl2805061 

116553025 

APOAl 

0.2718 


X 


X 




X 


rs651821 

116662579 

APOAl 

0.0845 


X 

X 

X 

X 

X 

X 



rs64 67 7 6 

109818530 

CELSR2 

0.2171 






X 




rs34216426 

56913088 

CETP 

0.0002 


X 


X 

X 


X 



rs5801 

56913513 

CETP 

0.1818 


X 


X 

X 

X 

X 



rs5802 

56919235 

CETP 

0.1194 


X 








rs3764261 

56993324 

CETP 

0.2719 

X 

X 

X 

X 

X 

X 

X 

X 


rs5883 

57007353 

CETP 

0.0401 

X 


X 

X 

X 

X 

X 

X 


rs5880 

57015091 

CETP 

0.0242 

X 

X 

X 

X 

X 

X 

X 

X 


rs2303790 

57017292 

CETP 

0.0006 

X 



X 

X 

X 

X 

X 

v_ 

_cl6_57095439 

57095439 

CETP 

0.0002 


X 


X 






rs509360 

61548559 

EADSl 

0.3364 








X 


rsl74546 

61569830 

EADSl 

0.4212 






X 




rs611229 

230324067 

GALNT2 

0.4077 


X 

X 

X 

X 

X 

X 

X 


rsl532085 

58683366 

LIPC 

0.4431 

X 

X 

X 

X 

X 

X 

X 

X 


rs261336 

58742418 

LIPC 

0.2302 

X 

X 

X 

X 

X 

X 

X 



rs28933094 

58855748 

LIPC 

0.0150 

X 

X 

X 

X 

X 

X 

X 

X 


rs268 

19813529 

LPL 

0.0180 

X 

X 

X 

X 

X 

X 

X 

X 


rs328 

19819724 

LPL 

0.0882 

X 

X 

X 

X 

X 

X 

X 

X 


rsl0096633 

19830921 

LPL 

0.0943 








X 

v_ 

_cll_47400044 

47400044 

MADD 

0.0005 

X 

X 

X 

X 

X 


X 



rs7946766 

48004369 

MADD 

0.1965 


X 

X 

X 

X 

X 

X 

X 


rsl2314392 

110010866 

MVK 

0.4033 







X 

X 


Table 6: Comparison of selection results for LDL. 






BH 

BH 



Across 

Across 


variant 

position 

locus 

MAE 

full 

marg Lasso 

SLOPE 

basic 

Trait 

Sites 

PLOS 

rs6756629 

44065090 

ABCG8 

0.0877 

X 

X 

X 

X 

X 

X 

X 

rsl45756111 

44102302 

ABCG8 

0.0115 


X 





X 

rs651821 

116662579 

APOAl 

0.0845 


X 

X 


X 

X 

X 

rsll216267 

116952392 

APOAl 

0.4578 


X 

X 

X 

X 

X 

X 

rs64 67 7 6 

109818530 

CELSR2 

0.2171 

X 

X 

X 

X 

X 

X 

X 

rsl2445698 

56928216 

CETP 

0.1887 







X 

rs3764261 

56993324 

CETP 

0.2719 





X 



rsl74546 

61569830 

EADSl 

0.4212 

X 

X 

X 

X 

X 

X 

X 

rs611229 

230324067 

GALNT2 

0.4077 





X 



rsl2610185 

19721722 

NCAN 

0.0617 







X 

rs2304130 

19789528 

NCAN 

0.0610 





X 
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Table 7: Comparison of selection results for TG. 






BH 

BH 




Across 

Across 


variant 

position 

locus 

MAF 

full 

raarg 

Lasso 

SLOPE 

basic 

Trait 

Sites 

PLOS 

rsl2805061 

116553025 

APOAl 

0.2718 


X 






X 

rs2266788 

116660686 

APOAl 

0.0891 








X 

rs3135506 

116662407 

APOAl 

0.0597 

X 

X 

X 

X 

X 

X 

X 

X 

rs651821 

116662579 

APOAl 

0.0845 

X 

X 

X 

X 

X 

X 

X 


rsll216267 

116952392 

APOAl 

0.4578 


X 







rsl561140 

56864398 

CETP 

0.4762 








X 

rs6591657 

61434532 

FADSl 

0.1506 








X 

rsl74546 

61569830 

FADSl 

0.4212 






X 



rs4846930 

230346829 

GALNT2 

0.4050 








X 

rsl260326 

27730940 

GCKR 

0.3566 

X 

X 

X 

X 

X 

X 

X 

X 

rs261336 

58742418 

LIPC 

0.2302 


X 


X 

X 

X 

X 

X 

rs28933094 

58855748 

LIPC 

0.0150 

X 

X 


X 

X 

X 

X 

X 

rs268 

19813529 

LPL 

0.0180 


X 


X 

X 

X 

X 

X 

rs328 

19819724 

LPL 

0.0882 

X 

X 

X 

X 

X 

X 

X 

X 

rsl0096633 

19830921 

LPL 

0.0943 








X 

rs2304130 

19789528 

NCAN 

0.0610 


X 


X 

X 

X 

X 

X 


when only the Across Traits approach selected it for HDL. Again, this locus has previous evidence 
of association with TG. 

To conclude, we note that applying BH to each phenotype separately gives different results 
in some cases. When applied to the p-values from the full model for HDL only, BH selects four 
variants that it does not select when applied to all traits simultaneously: rs2066715 in ABGAl, 
rs651821 in APOAl, rs34216426 in GETP, and rs611229 in GALNT2. All of these are also selected 
for HDL by most other methods. When applied to marginal p-values for HDL only, BH selects 
three variants that it does not select when applied to all traits simultaneously: rs62136410 in 
ABGG8, rs2303790 in GETP, and rsll988 in MADD. Most methods select rs2303790 for HDL, but 
no other methods select the other two of these; and ABGG8 does not even have prior evidence 
of association to HDL. The final difference arising from applying BH to traits separately is that it 
does not select rsl45756111 in ABGG8 for association with LDL when using marginal p-values, 
which brings it into agreement with the other methods on this variant. 

Effect of pruning variants 

Eor the case study, we pruned the set of variants so that the maximum of the absolute values 
of the pairwise correlations is 0.3 in order to facilitate comparison of the selection results by 
different methods. We will now describe some experiments we performed with the actual data 
to investigate the effects of this pruning. Eirst, we compare our selection results presented in the 
main text and above to the results when using X obtained by setting Cmax in Eigure equal to 
0.5, 0.7 and 0.9. Eor clarity, let X^^^ and X^^i denote the different versions of X. After 

5,000,000 MGMG iterations for each of four chains, over 99% of AZ^t and of AW^, or AGg are less 
than 0.05 for all levels of pruning, so convergence still seems acceptable. 

The values of Zyt usually do not change very much for different versions of X; and when 
they do, they almost always decrease. More precisely, much less than 1% of the values of Z-ui 
increase by more than 0.05. Eor given prior and trait, the percentage of values that change by more 
than 0.05 is less than 3% between X^^i and X^^i and 1% between X^^i and X^^i and between X^^^ 
and Xi^i. 

Table [^compares the total number of variants selected by our Bayesian priors for each of the 
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Table 8: Comparison of selection results from BH and Bayesian methods applied to different 
levels of variant pruning. The columns labeled R give the number of variants selected, and 
the columns labeled BFDR report the Bayesian FDR. The selection threshold for the Bayesian 
methods was ^ = 0.7. For BH, a = 0.05. 



HDL 

LDL 

TG 

prior 

R BFDR 

R BFDR 

R BFDR 


maximum correlation 0.3 


basic 

21 0.046 

4 0.012 

8 0.011 

Across Traits 

19 0.084 

8 0.059 

9 0.029 

Across Sites 

25 0.064 

5 0.063 

8 0.007 

BH full 

13 

3 

5 

BH marginal 

22 

6 

10 


maximum correlation 0.5 


basic 

21 

0.068 

4 

0.014 

9 

0.024 

Across Traits 

18 

0.091 

8 

0.060 

10 

0.044 

Across Sites 

25 

0.076 

4 

0.012 

9 

0.020 

maximum correlation 0.7 

basic 

20 

0.067 

4 

0.014 

9 

0.033 

Across Traits 

17 

0.081 

8 

0.057 

9 

0.022 

Across Sites 

24 

0.071 

4 

0.013 

9 

0.029 

maximum correlation 0.9 

basic 

19 

0.073 

4 

0.038 

9 

0.050 

Across Traits 

15 

0.071 

8 

0.090 

9 

0.048 

Across Sites 

22 

0.071 

4 

0.037 

9 

0.045 

BH full 

6 


2 


4 


BH marginal 

32 


13 


15 



pruning levels; they still did not select any variants that are in a locus (CRY2, G6PC2, MTNRIB, 
and PANKl) lacking any prior evidence of association to lipid traits. The table also shows the 
number of variants selected by applying BH to the p-values from and For a more 

detailed look at the differences. Table lists all cases in which a variant would be selected for 
one version of X but not another. These differences are explained in detail at the end of this 
section, but for now we categorize them as follows: (1) one variant {rs2066715) is replaced in 
the selections by another strongly correlated variant {rs2853579) when it enters X^^); (2) two rare 
variants are not in X^^) but are selected (or close to the threshold) for the other versions of X, at 
least in part because one of the 2-3 subjects with the minor allele has an extreme value for the 
trait; (3) two variants whose Z^,t is barely above the selection threshold for some trait and prior 
with X^^) but then decreases by as much as 0.15 as correlated variants enter X with Z^f < 0.1; and 
(4) three variants that are selected for X^^) but their Z^t decreases by as much as 0.5 as correlated 
variants enter X with Z^t > 0.15. 

To summarize, the different versions of X considered here do not have much effect on the 
selections made when computing Z^t from MCMC across all loci simultaneously. Very few of the 
variants that were left out of X^^^ are selected when they are included in a different version of X. 
In fact, adding correlated variants to X is more likely to result in a slight reduction in the number 
of selected variants. In contrast, applying BH to the p-values has a much greater effect on the 
selections between X^^^ and X^^^, with marginal p-values resulting in many more selections and 
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Table 9: Cases in which a variant would be selected for one version of X but not another. Entries 
show Zi,f; "NA" means that the variant is not in that version of X, while blanks mean that the 
variant would not be selected for any X. 



none 

Across Traits 

Across Sites 

max correlation 

0.3 0.5 0.7 0.9 

0.3 0.5 0.7 0.9 

0.3 0.5 0.7 0.9 


HDL 


rs2066715 

0.92 

0.89 

0.23 

0.27 

0.78 

0.75 

0.17 

0.22 

0.92 

0.89 

0.13 

0.13 

rs2853579 

NA 

NA 

0.75 

0.70 

NA 

NA 

0.75 

0.69 

NA 

NA 

0.88 

0.87 

rsl40547417 









NA 

0.73 

0.74 

0.66 

rs2575875 

0.92 

0.73 

0.40 

0.40 

0.79 

0.66 

0.39 

0.41 

0.94 

0.82 

0.46 

0.47 

rsl2314392 









0.73 

0.54 

0.53 

0.19 

rs7946766 

0.99 

0.97 

0.97 

0.56 

0.98 

0.97 

0.97 

0.52 

1.00 

0.98 

0.98 

0.65 

rsl74546 





0.72 

0.70 

0.62 

0.58 





LDL 

rs651821 


always >0.99 

0.73 

0.66 

0.63 

0.64 


TG 


v_c9_107544165 

NA 0.86 0.85 0.83 

NA 

0.83 

0.82 

0.79 

NA 0.87 0.87 0.87 

rsl74546 


0.74 

0.74 

0.68 

0.64 



p-values from the full model resulting in many fewer selections. 

Variants with sufficiently strong correlations could, however, lead to convergence problems 
in our MCMC. To allow arbitrary correlations, we also compute the posterior distribution exactly 


instead of using MCMC. To achieve this, we make the following simplifications as in Servin and 
Stephens (20071; Hormozdiari et al. ( 2014[ |; Chen et ^ ( |2015| : 


• Specify t. We choose t = 0.05 because the mean of the sampled values of t in the MCMC 
was always between 0.04 and 0.06 for all versions of X with the actual phenotype data. 

• Assume the posterior density is zero if | Z | > K, where K is typically no more than six but 
may be as small as one. 


Even requiring |Z| <4, however, would still leave far too many subsets if we were to consider 
all 1302 variants simultaneously, so we only consider one locus at a time. In fact, we consider 
only two loci: ABCAl and CETP, both with trait HDL, because these show stronger evidence of 
influencing the trait via multiple variants. The locus ABCAl has 58 variants, so we compute the 
posterior for |Z| <6; whereas the locus CETP has 98 variants, so we require |Z| <5. Pinally, we 
only tried this with our basic model. 

In addition to estimating the marginal posterior expected values as E [Z^ |y] « 'E,z:Z„=i /z,t (Z, t 
0.05 |y), we also consider two other approaches to selection. One is simply to choose Z with the 


largest value of /z,t(Z, t = 0.05|y). The other is to use the confidence set proposed by Hor¬ 


mozdiari et al. ( |2014| |, described as follows. The probability that a set S of variants contains 
all causal variants is P(iS) = Ez:Z„= 0 fori; 65 /z,T(Z ,t = 0.05|y). Ideally, one would choose the 
smallest S such that 7^(5) is above a specified threshold, but this is not computationally feasible. 
Instead, we use the stepwise selection process in Hormozdiari et al. (|2014||. At each step, the 
variant that increases 'P{S) the most is added to S. (Hormozdiari et al. 2014| use the stopping 
criterion 'P{S) > 0.95, but this results in a very large number of selections for our data when the 
maximum value of |Z| is four or greater. Instead, we use threshold 0.7. 


42 



























































The 70%-confidence set contains all the variants selected by any of the Bayesian methods, so 
Table [10| lists the variants in the order in which they enter the 70%-confidence set. The selections 
from BH applied to the p-values for only the one locus and the HDL are also shown, although 
BH applied to the marginal p-values for locus CETP selects nine additional variants that are 
not shown. [Chen et ’aT. ( |2015| l give a good argument for preferring selections based on marginal 
posterior expectations over the confidence set if there is more than one causal variant, so we do 
not discuss the confidence sets further here. BH applied to marginal p-values is probably also 
selecting too many variants. The other methods highlight many of the same variants. Overall, 
for 90% of the variants in the value of computed by MCMC is within 0.01 of the exact 
posterior marginals without pruning, and the latter value is almost always less than the former, 
which again suggests that pruning did not result in many missed discoveries. One of the largest 
differences between these two values is for rs2853579, which—as noted above in connection with 
Table |9] —is strongly correlated with rs2066715. Furthermore, for locus ABCAl, the set Z of 
variants with the largest posterior density is equal to the variants selected by using with 
For the locus CFTP, these two sets cannot be the same because Z-^t selects six variants, while we 
had to assume |Z| < 5 in order to compute the exact distribution; additionally, the mode of the 
posterior density includes one very rare variant (MAP 0.0002) that is in X^^^ but that the basic 
prior does not select. 

For completeness, the detailed descriptions of the differences in Table are as follows: 


rs2066715 and rs2853579, which are both in ABCAl, are in the table for all three priors applied 
to response HDF. Their correlation is just below 0.7, and Z^t for the former decreases dras¬ 
tically when the latter is in X and is selected instead—although Z^t for rs2853579 falls just 
below the selection threshold with Across Traits applied to X^^). The former is a missense 
variant predicted to be benign; whereas the latter is a synonymous variant. 

rsl40547417 is not in X^^), has Z^t barely above the selection threshold for X^^^ and X^^^ when 
Across Sites is applied to HDF, then drops slightly below the threshold for X^^). Its minor 
allele occurs in only three subjects, one of which has the 32nd largest value of HDF and 
another is in the 15th-percentile. Furthermore, this variant is in gene CFTP, and its group 
includes a variant with Zj^t ~ 0.8 and one with Z^t ~ 0.3. 


rs2575875 is strongly correlated with a new variant in X^^) (with Z^jt = 0.24) and with another 
new variant in X^^^—at which point all three variants have Z^t in 0.3-0.4. 

rsl2314392 has Z^t barely above the selection threshold when Across Sites is applied to X^^^. It 
has correlation 0.3-0.4 with three variants that are in X^^^. Two of these variants have 
Zyt < 0.05, while the other has Z^t = 0.2; this is easily enough for rsl2314392 to drop 
below the selection threshold. None of the new variants in X^^^ are strongly correlated with 
rsl2314392, but four of the new variants in X^^) have correlations 0.43-0.84 with rsl2314392. 
Three of these have Z^t < 0.1; but the one with correlation 0.75 has Z^t = 0.60, causing Z^t 
for rsl2314392 to decrease even more. 


rs7946766 has Z^t > 0.96 for HDF with all three priors until several variants strongly correlated 
with it are added to X^^^, including two that have Z^t ~ 0.3. 

v_c9_107544165 has minor allele that occurs in only two subjects, one of which has the 2nd 
largest value of HDF and the other is in the 14th-percentile. It is in the UTR 3' of ABCAl, 
which means it is in a group by itself in the Across Sites prior. 
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Table 10: Comparison of selection methods to see effect of using all variants. These tables list all variants in the 70%-confidence set in 
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rs9282537 



rsl74546 and rs651821 both have Z^t between 0.7 and 0.75 for but it decreases as more 
variables are added to X until it is between 0.58 and 0.65 for X^^); in other words, these 
variants are near the threshold for selection for all versions of X. 
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